Assignment 10: Surface currents from the OSCAR model¶

Let's access and plot global surface currents from OSCAR: https://www.esr.org/research/oscar/oscar-surface-currents/ according to Earth and Space Research, who create this data product:

"The Ocean Surface Current Analyses Real-time (OSCAR) is a NASA funded research project and is a direct computation of global surface currents using satellite sea surface height, wind, and temperature. Currents are calculated using a quasi-steady geostrophic model together with an eddy viscosity based wind-driven ageostrophic component and a thermal wind adjustment. The model calculates a surface current averaged over the top 30 m of the upper ocean.

In [1]:
# Our usual imports
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import cool_maps.plot as cplt
import cartopy.crs as ccrs
import cmocean.cm as cmo
import numpy as np
import cartopy.feature as cfeature
%matplotlib inline

The data is available from the NASA Physical Oceanography Distributed Active Archive Center: https://podaac.jpl.nasa.gov/

We can find more information about this data product at https://podaac.jpl.nasa.gov/dataset/OSCAR_L4_OC_third-deg?ids=&values=&search=oscar&provider=PODAAC

image.png

Download the 2018 dataset off Canvas here. It used to be accessible at https://podaac-opendap.jpl.nasa.gov/opendap/allData/oscar/preview/L4/oscar_third_deg/oscar_vel2018.nc.gz

In [2]:
# Download this file off Canvas
local = '/Users/csmit/Downloads/oscar_vel2018.nc.gz'
ds  = xr.open_dataset(local)

ds
Out[2]:
<xarray.Dataset> Size: 1GB
Dimensions:    (time: 72, year: 72, depth: 1, latitude: 481, longitude: 1201)
Coordinates:
  * time       (time) datetime64[ns] 576B 2018-01-01 2018-01-06 ... 2018-12-26
  * year       (year) float32 288B 2.018e+03 2.018e+03 ... 2.019e+03 2.019e+03
  * depth      (depth) float32 4B 15.0
  * latitude   (latitude) float64 4kB 80.0 79.67 79.33 ... -79.33 -79.67 -80.0
  * longitude  (longitude) float64 10kB 20.0 20.33 20.67 ... 419.3 419.7 420.0
Data variables:
    u          (time, depth, latitude, longitude) float64 333MB ...
    v          (time, depth, latitude, longitude) float64 333MB ...
    um         (time, depth, latitude, longitude) float64 333MB ...
    vm         (time, depth, latitude, longitude) float64 333MB ...
Attributes: (12/15)
    VARIABLE:       Ocean Surface Currents
    DATATYPE:       1/72 YEAR Interval
    DATASUBTYPE:    unfiltered
    GEORANGE:       20 to 420 -80 to 80
    PERIOD:         Jan.01,2018 to Dec.26,2018
    year:           2018
    ...             ...
    source:         Gary Lagerloef, ESR (lager@esr.org) and Kathleen Dohan, E...
    contact:        Kathleen Dohan (kdohan@esr.org) or John T. Gunn (gunn@esr...
    company:        Earth & Space Research, Seattle, WA
    reference:      Bonjean F. and G.S.E. Lagerloef, 2002 ,Diagnostic model a...
    note1:          Maximum Mask velocity is the geostrophic component at all...
    note2:          Longitude extends from 20 E to 420 E to avoid a break in ...
xarray.Dataset
    • time: 72
    • year: 72
    • depth: 1
    • latitude: 481
    • longitude: 1201
    • time
      (time)
      datetime64[ns]
      2018-01-01 ... 2018-12-26
      long_name :
      Day since 1992-10-05 00:00:00
      array(['2018-01-01T00:00:00.000000000', '2018-01-06T00:00:00.000000000',
             '2018-01-11T00:00:00.000000000', '2018-01-16T00:00:00.000000000',
             '2018-01-21T00:00:00.000000000', '2018-01-26T00:00:00.000000000',
             '2018-01-31T00:00:00.000000000', '2018-02-05T00:00:00.000000000',
             '2018-02-10T00:00:00.000000000', '2018-02-15T00:00:00.000000000',
             '2018-02-20T00:00:00.000000000', '2018-02-25T00:00:00.000000000',
             '2018-03-02T00:00:00.000000000', '2018-03-07T00:00:00.000000000',
             '2018-03-12T00:00:00.000000000', '2018-03-18T00:00:00.000000000',
             '2018-03-23T00:00:00.000000000', '2018-03-28T00:00:00.000000000',
             '2018-04-02T00:00:00.000000000', '2018-04-07T00:00:00.000000000',
             '2018-04-12T00:00:00.000000000', '2018-04-17T00:00:00.000000000',
             '2018-04-22T00:00:00.000000000', '2018-04-27T00:00:00.000000000',
             '2018-05-02T00:00:00.000000000', '2018-05-07T00:00:00.000000000',
             '2018-05-12T00:00:00.000000000', '2018-05-17T00:00:00.000000000',
             '2018-05-22T00:00:00.000000000', '2018-05-28T00:00:00.000000000',
             '2018-06-02T00:00:00.000000000', '2018-06-07T00:00:00.000000000',
             '2018-06-12T00:00:00.000000000', '2018-06-17T00:00:00.000000000',
             '2018-06-22T00:00:00.000000000', '2018-06-27T00:00:00.000000000',
             '2018-07-02T00:00:00.000000000', '2018-07-07T00:00:00.000000000',
             '2018-07-12T00:00:00.000000000', '2018-07-17T00:00:00.000000000',
             '2018-07-22T00:00:00.000000000', '2018-07-27T00:00:00.000000000',
             '2018-08-01T00:00:00.000000000', '2018-08-06T00:00:00.000000000',
             '2018-08-12T00:00:00.000000000', '2018-08-17T00:00:00.000000000',
             '2018-08-22T00:00:00.000000000', '2018-08-27T00:00:00.000000000',
             '2018-09-01T00:00:00.000000000', '2018-09-06T00:00:00.000000000',
             '2018-09-11T00:00:00.000000000', '2018-09-16T00:00:00.000000000',
             '2018-09-21T00:00:00.000000000', '2018-09-26T00:00:00.000000000',
             '2018-10-01T00:00:00.000000000', '2018-10-06T00:00:00.000000000',
             '2018-10-11T00:00:00.000000000', '2018-10-16T00:00:00.000000000',
             '2018-10-22T00:00:00.000000000', '2018-10-27T00:00:00.000000000',
             '2018-11-01T00:00:00.000000000', '2018-11-06T00:00:00.000000000',
             '2018-11-11T00:00:00.000000000', '2018-11-16T00:00:00.000000000',
             '2018-11-21T00:00:00.000000000', '2018-11-26T00:00:00.000000000',
             '2018-12-01T00:00:00.000000000', '2018-12-06T00:00:00.000000000',
             '2018-12-11T00:00:00.000000000', '2018-12-16T00:00:00.000000000',
             '2018-12-21T00:00:00.000000000', '2018-12-26T00:00:00.000000000'],
            dtype='datetime64[ns]')
    • year
      (year)
      float32
      2.018e+03 2.018e+03 ... 2.019e+03
      units :
      time in years
      long_name :
      Time in fractional year
      array([2018.    , 2018.0139, 2018.0278, 2018.0416, 2018.0555, 2018.0695,
             2018.0834, 2018.0972, 2018.1111, 2018.125 , 2018.1389, 2018.1528,
             2018.1666, 2018.1805, 2018.1945, 2018.2084, 2018.2222, 2018.2361,
             2018.25  , 2018.2639, 2018.2778, 2018.2916, 2018.3055, 2018.3195,
             2018.3334, 2018.3472, 2018.3611, 2018.375 , 2018.3889, 2018.4028,
             2018.4166, 2018.4305, 2018.4445, 2018.4584, 2018.4722, 2018.4861,
             2018.5   , 2018.5139, 2018.5278, 2018.5416, 2018.5555, 2018.5695,
             2018.5834, 2018.5972, 2018.6111, 2018.625 , 2018.6389, 2018.6528,
             2018.6666, 2018.6805, 2018.6945, 2018.7084, 2018.7222, 2018.7361,
             2018.75  , 2018.7639, 2018.7778, 2018.7916, 2018.8055, 2018.8195,
             2018.8334, 2018.8472, 2018.8611, 2018.875 , 2018.8889, 2018.9028,
             2018.9166, 2018.9305, 2018.9445, 2018.9584, 2018.9722, 2018.9861],
            dtype=float32)
    • depth
      (depth)
      float32
      15.0
      units :
      meter
      long_name :
      Depth
      array([15.], dtype=float32)
    • latitude
      (latitude)
      float64
      80.0 79.67 79.33 ... -79.67 -80.0
      units :
      degrees-north
      long_name :
      Latitude
      array([ 80.      ,  79.666667,  79.333333, ..., -79.333333, -79.666667,
             -80.      ])
    • longitude
      (longitude)
      float64
      20.0 20.33 20.67 ... 419.7 420.0
      units :
      degrees-east
      long_name :
      Longitude
      array([ 20.      ,  20.333333,  20.666667, ..., 419.333333, 419.666667,
             420.      ])
    • u
      (time, depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Zonal Currents
      [41593032 values with dtype=float64]
    • v
      (time, depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Meridional Currents
      [41593032 values with dtype=float64]
    • um
      (time, depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Zonal Currents Maximum Mask
      [41593032 values with dtype=float64]
    • vm
      (time, depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Meridional Currents Maximum Mask
      [41593032 values with dtype=float64]
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-01-01', '2018-01-06', '2018-01-11', '2018-01-16',
                     '2018-01-21', '2018-01-26', '2018-01-31', '2018-02-05',
                     '2018-02-10', '2018-02-15', '2018-02-20', '2018-02-25',
                     '2018-03-02', '2018-03-07', '2018-03-12', '2018-03-18',
                     '2018-03-23', '2018-03-28', '2018-04-02', '2018-04-07',
                     '2018-04-12', '2018-04-17', '2018-04-22', '2018-04-27',
                     '2018-05-02', '2018-05-07', '2018-05-12', '2018-05-17',
                     '2018-05-22', '2018-05-28', '2018-06-02', '2018-06-07',
                     '2018-06-12', '2018-06-17', '2018-06-22', '2018-06-27',
                     '2018-07-02', '2018-07-07', '2018-07-12', '2018-07-17',
                     '2018-07-22', '2018-07-27', '2018-08-01', '2018-08-06',
                     '2018-08-12', '2018-08-17', '2018-08-22', '2018-08-27',
                     '2018-09-01', '2018-09-06', '2018-09-11', '2018-09-16',
                     '2018-09-21', '2018-09-26', '2018-10-01', '2018-10-06',
                     '2018-10-11', '2018-10-16', '2018-10-22', '2018-10-27',
                     '2018-11-01', '2018-11-06', '2018-11-11', '2018-11-16',
                     '2018-11-21', '2018-11-26', '2018-12-01', '2018-12-06',
                     '2018-12-11', '2018-12-16', '2018-12-21', '2018-12-26'],
                    dtype='datetime64[ns]', name='time', freq=None))
    • year
      PandasIndex
      PandasIndex(Index([            2018.0,  2018.013916015625,   2018.02783203125,
             2018.0416259765625, 2018.0555419921875, 2018.0694580078125,
             2018.0833740234375,   2018.09716796875,  2018.111083984375,
                       2018.125,  2018.138916015625,   2018.15283203125,
             2018.1666259765625, 2018.1805419921875, 2018.1944580078125,
             2018.2083740234375,   2018.22216796875,  2018.236083984375,
                        2018.25,  2018.263916015625,   2018.27783203125,
             2018.2916259765625, 2018.3055419921875, 2018.3194580078125,
             2018.3333740234375,   2018.34716796875,  2018.361083984375,
                       2018.375,  2018.388916015625,   2018.40283203125,
             2018.4166259765625, 2018.4305419921875, 2018.4444580078125,
             2018.4583740234375,   2018.47216796875,  2018.486083984375,
                         2018.5,  2018.513916015625,   2018.52783203125,
             2018.5416259765625, 2018.5555419921875, 2018.5694580078125,
             2018.5833740234375,   2018.59716796875,  2018.611083984375,
                       2018.625,  2018.638916015625,   2018.65283203125,
             2018.6666259765625, 2018.6805419921875, 2018.6944580078125,
             2018.7083740234375,   2018.72216796875,  2018.736083984375,
                        2018.75,  2018.763916015625,   2018.77783203125,
             2018.7916259765625, 2018.8055419921875, 2018.8194580078125,
             2018.8333740234375,   2018.84716796875,  2018.861083984375,
                       2018.875,  2018.888916015625,   2018.90283203125,
             2018.9166259765625, 2018.9305419921875, 2018.9444580078125,
             2018.9583740234375,   2018.97216796875,  2018.986083984375],
            dtype='float32', name='year'))
    • depth
      PandasIndex
      PandasIndex(Index([15.0], dtype='float32', name='depth'))
    • latitude
      PandasIndex
      PandasIndex(Index([              80.0,  79.66666666666667,  79.33333333333333,
                           79.0,  78.66666666666667,  78.33333333333333,
                           78.0,  77.66666666666667,  77.33333333333333,
                           77.0,
             ...
                          -77.0, -77.33333333333333, -77.66666666666667,
                          -78.0, -78.33333333333333, -78.66666666666667,
                          -79.0, -79.33333333333333, -79.66666666666667,
                          -80.0],
            dtype='float64', name='latitude', length=481))
    • longitude
      PandasIndex
      PandasIndex(Index([              20.0, 20.333333333333332, 20.666666666666664,
                           21.0, 21.333333333333332, 21.666666666666664,
                           22.0, 22.333333333333332, 22.666666666666664,
                           23.0,
             ...
                          417.0,  417.3333333333333,  417.6666666666667,
                          418.0,  418.3333333333333,  418.6666666666667,
                          419.0,  419.3333333333333,  419.6666666666667,
                          420.0],
            dtype='float64', name='longitude', length=1201))
  • VARIABLE :
    Ocean Surface Currents
    DATATYPE :
    1/72 YEAR Interval
    DATASUBTYPE :
    unfiltered
    GEORANGE :
    20 to 420 -80 to 80
    PERIOD :
    Jan.01,2018 to Dec.26,2018
    year :
    2018
    description :
    OSCAR Third Degree Sea Surface Velocity
    CREATION_DATE :
    03:39 30-Jan-2019
    version :
    2009.0
    source :
    Gary Lagerloef, ESR (lager@esr.org) and Kathleen Dohan, ESR (kdohan@esr.org)
    contact :
    Kathleen Dohan (kdohan@esr.org) or John T. Gunn (gunn@esr.org)
    company :
    Earth & Space Research, Seattle, WA
    reference :
    Bonjean F. and G.S.E. Lagerloef, 2002 ,Diagnostic model and analysis of the surface currents in the tropical Pacific ocean, J. Phys. Oceanogr., 32, 2,938-2,954
    note1 :
    Maximum Mask velocity is the geostrophic component at all points + any concurrent Ekman and buoyancy components
    note2 :
    Longitude extends from 20 E to 420 E to avoid a break in major ocean basins. Data repeats in overlap region.

Let's first plot global surface U (zonal speed) on a particular time record: 2018-01-01¶

In [3]:
ds.u.sel(time='2018-01-01', method='nearest').plot();
No description has been provided for this image

So we can see it is a global dataset. How about we just look at the western north atlantic?¶

Select only lon = 275:350, lat = 60:20, for a single time record: 2018-01-01 remember to use .sel( dim = slice(start, stop)

In [4]:
ds.u.sel(time='2018-01-01', longitude =slice(275,350), latitude=slice(60,20) ).plot( vmin=-.75 , vmax=.75, cmap='RdBu');
No description has been provided for this image
In [5]:
ds.v.sel(time='2018-01-01', longitude =slice(275,350), latitude=slice(60,20) ).plot( vmin=-.75 , vmax=.75, cmap='RdBu');
No description has been provided for this image

Let's save that exact subset you just created into a new dataarray called subset. This will streamline the coding as we move forward, rather than calling the entire dataset each time.¶

In [6]:
subset = ds.sel(time='2018-01-01', longitude =slice(275,350), latitude=slice(60,20) )
subset
Out[6]:
<xarray.Dataset> Size: 878kB
Dimensions:    (year: 72, depth: 1, latitude: 121, longitude: 226)
Coordinates:
    time       datetime64[ns] 8B 2018-01-01
  * year       (year) float32 288B 2.018e+03 2.018e+03 ... 2.019e+03 2.019e+03
  * depth      (depth) float32 4B 15.0
  * latitude   (latitude) float64 968B 60.0 59.67 59.33 ... 20.67 20.33 20.0
  * longitude  (longitude) float64 2kB 275.0 275.3 275.7 ... 349.3 349.7 350.0
Data variables:
    u          (depth, latitude, longitude) float64 219kB ...
    v          (depth, latitude, longitude) float64 219kB ...
    um         (depth, latitude, longitude) float64 219kB ...
    vm         (depth, latitude, longitude) float64 219kB ...
Attributes: (12/15)
    VARIABLE:       Ocean Surface Currents
    DATATYPE:       1/72 YEAR Interval
    DATASUBTYPE:    unfiltered
    GEORANGE:       20 to 420 -80 to 80
    PERIOD:         Jan.01,2018 to Dec.26,2018
    year:           2018
    ...             ...
    source:         Gary Lagerloef, ESR (lager@esr.org) and Kathleen Dohan, E...
    contact:        Kathleen Dohan (kdohan@esr.org) or John T. Gunn (gunn@esr...
    company:        Earth & Space Research, Seattle, WA
    reference:      Bonjean F. and G.S.E. Lagerloef, 2002 ,Diagnostic model a...
    note1:          Maximum Mask velocity is the geostrophic component at all...
    note2:          Longitude extends from 20 E to 420 E to avoid a break in ...
xarray.Dataset
    • year: 72
    • depth: 1
    • latitude: 121
    • longitude: 226
    • time
      ()
      datetime64[ns]
      2018-01-01
      long_name :
      Day since 1992-10-05 00:00:00
      array('2018-01-01T00:00:00.000000000', dtype='datetime64[ns]')
    • year
      (year)
      float32
      2.018e+03 2.018e+03 ... 2.019e+03
      units :
      time in years
      long_name :
      Time in fractional year
      array([2018.    , 2018.0139, 2018.0278, 2018.0416, 2018.0555, 2018.0695,
             2018.0834, 2018.0972, 2018.1111, 2018.125 , 2018.1389, 2018.1528,
             2018.1666, 2018.1805, 2018.1945, 2018.2084, 2018.2222, 2018.2361,
             2018.25  , 2018.2639, 2018.2778, 2018.2916, 2018.3055, 2018.3195,
             2018.3334, 2018.3472, 2018.3611, 2018.375 , 2018.3889, 2018.4028,
             2018.4166, 2018.4305, 2018.4445, 2018.4584, 2018.4722, 2018.4861,
             2018.5   , 2018.5139, 2018.5278, 2018.5416, 2018.5555, 2018.5695,
             2018.5834, 2018.5972, 2018.6111, 2018.625 , 2018.6389, 2018.6528,
             2018.6666, 2018.6805, 2018.6945, 2018.7084, 2018.7222, 2018.7361,
             2018.75  , 2018.7639, 2018.7778, 2018.7916, 2018.8055, 2018.8195,
             2018.8334, 2018.8472, 2018.8611, 2018.875 , 2018.8889, 2018.9028,
             2018.9166, 2018.9305, 2018.9445, 2018.9584, 2018.9722, 2018.9861],
            dtype=float32)
    • depth
      (depth)
      float32
      15.0
      units :
      meter
      long_name :
      Depth
      array([15.], dtype=float32)
    • latitude
      (latitude)
      float64
      60.0 59.67 59.33 ... 20.33 20.0
      units :
      degrees-north
      long_name :
      Latitude
      array([60.      , 59.666667, 59.333333, 59.      , 58.666667, 58.333333,
             58.      , 57.666667, 57.333333, 57.      , 56.666667, 56.333333,
             56.      , 55.666667, 55.333333, 55.      , 54.666667, 54.333333,
             54.      , 53.666667, 53.333333, 53.      , 52.666667, 52.333333,
             52.      , 51.666667, 51.333333, 51.      , 50.666667, 50.333333,
             50.      , 49.666667, 49.333333, 49.      , 48.666667, 48.333333,
             48.      , 47.666667, 47.333333, 47.      , 46.666667, 46.333333,
             46.      , 45.666667, 45.333333, 45.      , 44.666667, 44.333333,
             44.      , 43.666667, 43.333333, 43.      , 42.666667, 42.333333,
             42.      , 41.666667, 41.333333, 41.      , 40.666667, 40.333333,
             40.      , 39.666667, 39.333333, 39.      , 38.666667, 38.333333,
             38.      , 37.666667, 37.333333, 37.      , 36.666667, 36.333333,
             36.      , 35.666667, 35.333333, 35.      , 34.666667, 34.333333,
             34.      , 33.666667, 33.333333, 33.      , 32.666667, 32.333333,
             32.      , 31.666667, 31.333333, 31.      , 30.666667, 30.333333,
             30.      , 29.666667, 29.333333, 29.      , 28.666667, 28.333333,
             28.      , 27.666667, 27.333333, 27.      , 26.666667, 26.333333,
             26.      , 25.666667, 25.333333, 25.      , 24.666667, 24.333333,
             24.      , 23.666667, 23.333333, 23.      , 22.666667, 22.333333,
             22.      , 21.666667, 21.333333, 21.      , 20.666667, 20.333333,
             20.      ])
    • longitude
      (longitude)
      float64
      275.0 275.3 275.7 ... 349.7 350.0
      units :
      degrees-east
      long_name :
      Longitude
      array([275.      , 275.333333, 275.666667, ..., 349.333333, 349.666667,
             350.      ])
    • u
      (depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Zonal Currents
      [27346 values with dtype=float64]
    • v
      (depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Meridional Currents
      [27346 values with dtype=float64]
    • um
      (depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Zonal Currents Maximum Mask
      [27346 values with dtype=float64]
    • vm
      (depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Meridional Currents Maximum Mask
      [27346 values with dtype=float64]
    • year
      PandasIndex
      PandasIndex(Index([            2018.0,  2018.013916015625,   2018.02783203125,
             2018.0416259765625, 2018.0555419921875, 2018.0694580078125,
             2018.0833740234375,   2018.09716796875,  2018.111083984375,
                       2018.125,  2018.138916015625,   2018.15283203125,
             2018.1666259765625, 2018.1805419921875, 2018.1944580078125,
             2018.2083740234375,   2018.22216796875,  2018.236083984375,
                        2018.25,  2018.263916015625,   2018.27783203125,
             2018.2916259765625, 2018.3055419921875, 2018.3194580078125,
             2018.3333740234375,   2018.34716796875,  2018.361083984375,
                       2018.375,  2018.388916015625,   2018.40283203125,
             2018.4166259765625, 2018.4305419921875, 2018.4444580078125,
             2018.4583740234375,   2018.47216796875,  2018.486083984375,
                         2018.5,  2018.513916015625,   2018.52783203125,
             2018.5416259765625, 2018.5555419921875, 2018.5694580078125,
             2018.5833740234375,   2018.59716796875,  2018.611083984375,
                       2018.625,  2018.638916015625,   2018.65283203125,
             2018.6666259765625, 2018.6805419921875, 2018.6944580078125,
             2018.7083740234375,   2018.72216796875,  2018.736083984375,
                        2018.75,  2018.763916015625,   2018.77783203125,
             2018.7916259765625, 2018.8055419921875, 2018.8194580078125,
             2018.8333740234375,   2018.84716796875,  2018.861083984375,
                       2018.875,  2018.888916015625,   2018.90283203125,
             2018.9166259765625, 2018.9305419921875, 2018.9444580078125,
             2018.9583740234375,   2018.97216796875,  2018.986083984375],
            dtype='float32', name='year'))
    • depth
      PandasIndex
      PandasIndex(Index([15.0], dtype='float32', name='depth'))
    • latitude
      PandasIndex
      PandasIndex(Index([              60.0,  59.66666666666667, 59.333333333333336,
                           59.0,  58.66666666666667, 58.333333333333336,
                           58.0,  57.66666666666667, 57.333333333333336,
                           57.0,
             ...
                           23.0,  22.66666666666667, 22.333333333333343,
                           22.0,  21.66666666666667, 21.333333333333343,
                           21.0,  20.66666666666667, 20.333333333333343,
                           20.0],
            dtype='float64', name='latitude', length=121))
    • longitude
      PandasIndex
      PandasIndex(Index([             275.0, 275.33333333333337,  275.6666666666667,
                          276.0, 276.33333333333337,  276.6666666666667,
                          277.0, 277.33333333333337,  277.6666666666667,
                          278.0,
             ...
                          347.0,  347.3333333333333,  347.6666666666667,
                          348.0,  348.3333333333333,  348.6666666666667,
                          349.0,  349.3333333333333,  349.6666666666667,
                          350.0],
            dtype='float64', name='longitude', length=226))
  • VARIABLE :
    Ocean Surface Currents
    DATATYPE :
    1/72 YEAR Interval
    DATASUBTYPE :
    unfiltered
    GEORANGE :
    20 to 420 -80 to 80
    PERIOD :
    Jan.01,2018 to Dec.26,2018
    year :
    2018
    description :
    OSCAR Third Degree Sea Surface Velocity
    CREATION_DATE :
    03:39 30-Jan-2019
    version :
    2009.0
    source :
    Gary Lagerloef, ESR (lager@esr.org) and Kathleen Dohan, ESR (kdohan@esr.org)
    contact :
    Kathleen Dohan (kdohan@esr.org) or John T. Gunn (gunn@esr.org)
    company :
    Earth & Space Research, Seattle, WA
    reference :
    Bonjean F. and G.S.E. Lagerloef, 2002 ,Diagnostic model and analysis of the surface currents in the tropical Pacific ocean, J. Phys. Oceanogr., 32, 2,938-2,954
    note1 :
    Maximum Mask velocity is the geostrophic component at all points + any concurrent Ekman and buoyancy components
    note2 :
    Longitude extends from 20 E to 420 E to avoid a break in major ocean basins. Data repeats in overlap region.

So we can visualize u and v, but how can we combine them? By plotting the vectors!¶

We can use the matplotlib function called .quiver() to plot vectors. It takes as arguments the x- and y- locations, and the u and v components of the velocity: plt.quiver(x, y, u, v)

The velocity components in our data have 3 dims: depth, latitude, and longitude, and quiver can only handle 2D data, so we need to select a dimension using the numpy style selections. In our case the u and v data we want for plotting can be gotten with subset.u[0,:,:] and subset.v[0,:,:]. This is useful notation if the third dimension is more than 1 count, since you can specify which index you want directly without having to re-subset the data.

In [7]:
plt.quiver( subset.longitude, subset.latitude, subset.u[0, :, :], subset.v[0,:,:] );
No description has been provided for this image

How about we zoom in, now that we've plotted the full (albeit subsetted) map, using plt.xlim() and plt.ylim()¶

In [8]:
plt.quiver( subset.longitude, subset.latitude, subset.u[0, :, :], subset.v[0,:,:] );
plt.xlim([280, 300]);
plt.ylim([25, 45]);
No description has been provided for this image

.quiver works just like .scatter in that you can specify an array for the colors as well¶

To do this, you need to get the magnitude of the velocity. If you picture a single u and v pair as sides of a right triangle: the u is your adjacent side, the v is your opposite side, how do we get the magnitude or the hypotenuse?¶

In [9]:
plt.quiver( subset.longitude, subset.latitude, subset.u[0, :, :], subset.v[0,:,:], (subset.u[0, :, :]**2 + subset.v[0, :, :]**2 )**0.5 )
plt.xlim([280, 300]);
plt.ylim([25, 45]);
plt.colorbar();
No description has been provided for this image

Now we know how our data looks, without any bells & whistles. This should always be your first step when working with a new dataset. Plot the data as simply as possible to make sure it looks right, before starting on more complicated visualizations.¶

We've done our data check, so we can move on and plot our data on top of a nice map from Mike Smith's cool_maps. Note cool_maps may struggle with the global visuals, since it's designed for more regional focus.¶

We have to do a bit more work with the latitude and longitude values to make this work.¶
In [10]:
print(subset.longitude.shape,subset.latitude.shape,subset.u.shape,subset.v.shape)
(226,) (121,) (1, 121, 226) (1, 121, 226)

The depth dimension in u and v is only 1 count, so it doesn't lend much to the structure. How do we get rid of extra dimensions again?¶

.squeeze¶

In [11]:
subset = subset.squeeze()
print(subset.longitude.shape,subset.latitude.shape,subset.u.shape,subset.v.shape)
(226,) (121,) (121, 226) (121, 226)

While having that mismatch of shapes works fine for a regular quiver plot because Python can figure out what needs to be matched to make it work, when we plot with cartopy or Mike Smith's cool_maps, we have to use transform. Doing so makes it too hard for Python to guess what needs to be done with the latitude and longitude to make them match the u and v. So we have to be a bit more detailed in our coding.¶

In [12]:
# First we turn lon and lat into variables with shape (121,226). 
# The longitude defaulted to 0:360, but cool_maps works with -180:180, so we subtract it by 360.
lon, lat = np.meshgrid(subset.longitude.values - 360,subset.latitude.values)

print('Before: ',lon.shape,lat.shape,subset.u.shape,subset.v.shape)

# Then we have to turn them all into 1D arrays for the transform to work properly. reshape(-1) does just that!
lon = lon.reshape(-1)
lat = lat.reshape(-1)
u   = subset.u.values.reshape(-1)
v   = subset.v.values.reshape(-1)
umag= (u**2 + v**2 )**0.5 

print('After: ',lon.shape,lat.shape,u.shape,v.shape)
Before:  (121, 226) (121, 226) (121, 226) (121, 226)
After:  (27346,) (27346,) (27346,) (27346,)
In [13]:
# Now we set the cool_maps extent
extent     = [-80, -40, 20, 45]

# Then borrow our cool_maps create line from past coding exercises
fig,ax  = cplt.create(extent, proj=ccrs.Mercator(), figsize=(10,6), oceancolor='w')
ax_posi = plt.quiver(lon,lat,u,v,umag,transform=ccrs.PlateCarree())

# Add in a title
plt.title('OSCAR Surface Currents on January 1, 2018');

# And a colorbar
plt.colorbar(label='Speed [m/s]');
No description has been provided for this image

Can we make a streamplot as well?¶

https://matplotlib.org/stable/plot_types/arrays/streamplot.html

We need to np.flip() our latitude, because it goes from 60 N to 20 N, but .streamplot needs values for x and y to be monotonically increasing. Since we np.flip() our latitude, we need to np.flipud() our u and v as well.

In [14]:
lon, lat = np.meshgrid(subset.longitude - 360,np.flip(subset.latitude))
u   = np.flipud(subset.u.values.squeeze())
v   = np.flipud(subset.v.values.squeeze())
umag= (u**2 + v**2 )**0.5 

plt.streamplot(lon,lat,u,v,color=umag,linewidth=umag,density=10,arrowsize=.5);
plt.colorbar();
No description has been provided for this image

#1. What specific products were used to calculate the OSCAR 1/3 degree product?¶

(More specific than the immediate description for the dataset)¶

hh

#2. Plot the annual mean of the surface currents (both global and regional subset) for 2018 as colored vectors.¶

Be sure to check the documentation on xarray means to ensure you take it along the correct dimension and skip NaNs.¶

additional imports¶

In [15]:
# Select surface level (assuming depth=15.0 is the only level)
surface_currents = ds.sel(depth=15.0)

# Calculate the annual mean for 2018
u_mean = surface_currents['u'].mean(dim='time', skipna=True)
v_mean = surface_currents['v'].mean(dim='time', skipna=True)

# Calculate magnitude of the currents
umag_g =(u_mean**2 + v_mean**2)**0.5

# Global plot
plt.figure(figsize=(14, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()  # Add coastlines
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')  # White background for land
ax.add_feature(cfeature.OCEAN, color='white')  # White background for oceans

quiver = plt.quiver(surface_currents['longitude'], surface_currents['latitude'], u_mean, v_mean, umag_g, transform=ccrs.PlateCarree(), scale=25, pivot='middle', cmap='viridis')

ax.set_xticks(np.arange(-180, 181, 30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90, 91, 30), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x}°"))
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda y, _: f"{y}°"))
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")

plt.colorbar(quiver, label='Surface Current Magnitude (m/s)', orientation='vertical')
plt.title('Annual Mean Surface Currents (Global) - 2018')
plt.show()

# Regional subset for the northwestern Atlantic
region = surface_currents.sel(latitude=slice(45, 25), longitude=slice(280, 300))
u_mean_region = region['u'].mean(dim='time', skipna=True)
v_mean_region = region['v'].mean(dim='time', skipna=True)
umag_r = (u_mean_region**2 + v_mean_region**2)**0.5

# Regional plot
plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([280, 300, 25, 45], crs=ccrs.PlateCarree())
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='white')
ax.set_xlabel('longitude')
quiver = plt.quiver(region['longitude'], region['latitude'], u_mean_region, v_mean_region, umag_r, transform=ccrs.PlateCarree(), scale=10, pivot='middle', cmap='viridis')

cbar = plt.colorbar(quiver, label='Surface Current Magnitude (m/s)', orientation='vertical')


# Add gridlines with proper longitude and latitude labels
gl = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}

plt.title('Annual Mean Surface Currents (Northwestern Atlantic) - 2018', fontsize=14)

# Adjust layout to prevent squishing
#plt.tight_layout()
plt.show()
No description has been provided for this image
No description has been provided for this image

#3. Plot the annual mean of the surface currents (both global and regional subset) for 2017 (available off Canvas as well) as colored vectors.¶

In [16]:
local = '/Users/csmit/Downloads/oscar_vel2017.nc.gz'
ds2  = xr.open_dataset(local)

ds2
Out[16]:
<xarray.Dataset> Size: 1GB
Dimensions:    (time: 72, year: 72, depth: 1, latitude: 481, longitude: 1201)
Coordinates:
  * time       (time) datetime64[ns] 576B 2017-01-01 2017-01-06 ... 2017-12-26
  * year       (year) float32 288B 2.017e+03 2.017e+03 ... 2.018e+03 2.018e+03
  * depth      (depth) float32 4B 15.0
  * latitude   (latitude) float64 4kB 80.0 79.67 79.33 ... -79.33 -79.67 -80.0
  * longitude  (longitude) float64 10kB 20.0 20.33 20.67 ... 419.3 419.7 420.0
Data variables:
    u          (time, depth, latitude, longitude) float64 333MB ...
    v          (time, depth, latitude, longitude) float64 333MB ...
    um         (time, depth, latitude, longitude) float64 333MB ...
    vm         (time, depth, latitude, longitude) float64 333MB ...
Attributes: (12/15)
    VARIABLE:       Ocean Surface Currents
    DATATYPE:       1/72 YEAR Interval
    DATASUBTYPE:    unfiltered
    GEORANGE:       20 to 420 -80 to 80
    PERIOD:         Jan.01,2017 to Dec.26,2017
    year:           2017
    ...             ...
    source:         Gary Lagerloef, ESR (lager@esr.org) and Kathleen Dohan, E...
    contact:        Kathleen Dohan (kdohan@esr.org) or John T. Gunn (gunn@esr...
    company:        Earth & Space Research, Seattle, WA
    reference:      Bonjean F. and G.S.E. Lagerloef, 2002 ,Diagnostic model a...
    note1:          Maximum Mask velocity is the geostrophic component at all...
    note2:          Longitude extends from 20 E to 420 E to avoid a break in ...
xarray.Dataset
    • time: 72
    • year: 72
    • depth: 1
    • latitude: 481
    • longitude: 1201
    • time
      (time)
      datetime64[ns]
      2017-01-01 ... 2017-12-26
      long_name :
      Day since 1992-10-05 00:00:00
      array(['2017-01-01T00:00:00.000000000', '2017-01-06T00:00:00.000000000',
             '2017-01-11T00:00:00.000000000', '2017-01-16T00:00:00.000000000',
             '2017-01-21T00:00:00.000000000', '2017-01-26T00:00:00.000000000',
             '2017-01-31T00:00:00.000000000', '2017-02-05T00:00:00.000000000',
             '2017-02-10T00:00:00.000000000', '2017-02-15T00:00:00.000000000',
             '2017-02-20T00:00:00.000000000', '2017-02-25T00:00:00.000000000',
             '2017-03-02T00:00:00.000000000', '2017-03-07T00:00:00.000000000',
             '2017-03-12T00:00:00.000000000', '2017-03-18T00:00:00.000000000',
             '2017-03-23T00:00:00.000000000', '2017-03-28T00:00:00.000000000',
             '2017-04-02T00:00:00.000000000', '2017-04-07T00:00:00.000000000',
             '2017-04-12T00:00:00.000000000', '2017-04-17T00:00:00.000000000',
             '2017-04-22T00:00:00.000000000', '2017-04-27T00:00:00.000000000',
             '2017-05-02T00:00:00.000000000', '2017-05-07T00:00:00.000000000',
             '2017-05-12T00:00:00.000000000', '2017-05-17T00:00:00.000000000',
             '2017-05-22T00:00:00.000000000', '2017-05-28T00:00:00.000000000',
             '2017-06-02T00:00:00.000000000', '2017-06-07T00:00:00.000000000',
             '2017-06-12T00:00:00.000000000', '2017-06-17T00:00:00.000000000',
             '2017-06-22T00:00:00.000000000', '2017-06-27T00:00:00.000000000',
             '2017-07-02T00:00:00.000000000', '2017-07-07T00:00:00.000000000',
             '2017-07-12T00:00:00.000000000', '2017-07-17T00:00:00.000000000',
             '2017-07-22T00:00:00.000000000', '2017-07-27T00:00:00.000000000',
             '2017-08-01T00:00:00.000000000', '2017-08-06T00:00:00.000000000',
             '2017-08-12T00:00:00.000000000', '2017-08-17T00:00:00.000000000',
             '2017-08-22T00:00:00.000000000', '2017-08-27T00:00:00.000000000',
             '2017-09-01T00:00:00.000000000', '2017-09-06T00:00:00.000000000',
             '2017-09-11T00:00:00.000000000', '2017-09-16T00:00:00.000000000',
             '2017-09-21T00:00:00.000000000', '2017-09-26T00:00:00.000000000',
             '2017-10-01T00:00:00.000000000', '2017-10-06T00:00:00.000000000',
             '2017-10-11T00:00:00.000000000', '2017-10-16T00:00:00.000000000',
             '2017-10-22T00:00:00.000000000', '2017-10-27T00:00:00.000000000',
             '2017-11-01T00:00:00.000000000', '2017-11-06T00:00:00.000000000',
             '2017-11-11T00:00:00.000000000', '2017-11-16T00:00:00.000000000',
             '2017-11-21T00:00:00.000000000', '2017-11-26T00:00:00.000000000',
             '2017-12-01T00:00:00.000000000', '2017-12-06T00:00:00.000000000',
             '2017-12-11T00:00:00.000000000', '2017-12-16T00:00:00.000000000',
             '2017-12-21T00:00:00.000000000', '2017-12-26T00:00:00.000000000'],
            dtype='datetime64[ns]')
    • year
      (year)
      float32
      2.017e+03 2.017e+03 ... 2.018e+03
      units :
      time in years
      long_name :
      Time in fractional year
      array([2017.    , 2017.0139, 2017.0278, 2017.0416, 2017.0555, 2017.0695,
             2017.0834, 2017.0972, 2017.1111, 2017.125 , 2017.1389, 2017.1528,
             2017.1666, 2017.1805, 2017.1945, 2017.2084, 2017.2222, 2017.2361,
             2017.25  , 2017.2639, 2017.2778, 2017.2916, 2017.3055, 2017.3195,
             2017.3334, 2017.3472, 2017.3611, 2017.375 , 2017.3889, 2017.4028,
             2017.4166, 2017.4305, 2017.4445, 2017.4584, 2017.4722, 2017.4861,
             2017.5   , 2017.5139, 2017.5278, 2017.5416, 2017.5555, 2017.5695,
             2017.5834, 2017.5972, 2017.6111, 2017.625 , 2017.6389, 2017.6528,
             2017.6666, 2017.6805, 2017.6945, 2017.7084, 2017.7222, 2017.7361,
             2017.75  , 2017.7639, 2017.7778, 2017.7916, 2017.8055, 2017.8195,
             2017.8334, 2017.8472, 2017.8611, 2017.875 , 2017.8889, 2017.9028,
             2017.9166, 2017.9305, 2017.9445, 2017.9584, 2017.9722, 2017.9861],
            dtype=float32)
    • depth
      (depth)
      float32
      15.0
      units :
      meter
      long_name :
      Depth
      array([15.], dtype=float32)
    • latitude
      (latitude)
      float64
      80.0 79.67 79.33 ... -79.67 -80.0
      units :
      degrees-north
      long_name :
      Latitude
      array([ 80.      ,  79.666667,  79.333333, ..., -79.333333, -79.666667,
             -80.      ])
    • longitude
      (longitude)
      float64
      20.0 20.33 20.67 ... 419.7 420.0
      units :
      degrees-east
      long_name :
      Longitude
      array([ 20.      ,  20.333333,  20.666667, ..., 419.333333, 419.666667,
             420.      ])
    • u
      (time, depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Zonal Currents
      [41593032 values with dtype=float64]
    • v
      (time, depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Meridional Currents
      [41593032 values with dtype=float64]
    • um
      (time, depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Zonal Currents Maximum Mask
      [41593032 values with dtype=float64]
    • vm
      (time, depth, latitude, longitude)
      float64
      ...
      units :
      meter/sec
      long_name :
      Ocean Surface Meridional Currents Maximum Mask
      [41593032 values with dtype=float64]
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2017-01-01', '2017-01-06', '2017-01-11', '2017-01-16',
                     '2017-01-21', '2017-01-26', '2017-01-31', '2017-02-05',
                     '2017-02-10', '2017-02-15', '2017-02-20', '2017-02-25',
                     '2017-03-02', '2017-03-07', '2017-03-12', '2017-03-18',
                     '2017-03-23', '2017-03-28', '2017-04-02', '2017-04-07',
                     '2017-04-12', '2017-04-17', '2017-04-22', '2017-04-27',
                     '2017-05-02', '2017-05-07', '2017-05-12', '2017-05-17',
                     '2017-05-22', '2017-05-28', '2017-06-02', '2017-06-07',
                     '2017-06-12', '2017-06-17', '2017-06-22', '2017-06-27',
                     '2017-07-02', '2017-07-07', '2017-07-12', '2017-07-17',
                     '2017-07-22', '2017-07-27', '2017-08-01', '2017-08-06',
                     '2017-08-12', '2017-08-17', '2017-08-22', '2017-08-27',
                     '2017-09-01', '2017-09-06', '2017-09-11', '2017-09-16',
                     '2017-09-21', '2017-09-26', '2017-10-01', '2017-10-06',
                     '2017-10-11', '2017-10-16', '2017-10-22', '2017-10-27',
                     '2017-11-01', '2017-11-06', '2017-11-11', '2017-11-16',
                     '2017-11-21', '2017-11-26', '2017-12-01', '2017-12-06',
                     '2017-12-11', '2017-12-16', '2017-12-21', '2017-12-26'],
                    dtype='datetime64[ns]', name='time', freq=None))
    • year
      PandasIndex
      PandasIndex(Index([            2017.0,  2017.013916015625,   2017.02783203125,
             2017.0416259765625, 2017.0555419921875, 2017.0694580078125,
             2017.0833740234375,   2017.09716796875,  2017.111083984375,
                       2017.125,  2017.138916015625,   2017.15283203125,
             2017.1666259765625, 2017.1805419921875, 2017.1944580078125,
             2017.2083740234375,   2017.22216796875,  2017.236083984375,
                        2017.25,  2017.263916015625,   2017.27783203125,
             2017.2916259765625, 2017.3055419921875, 2017.3194580078125,
             2017.3333740234375,   2017.34716796875,  2017.361083984375,
                       2017.375,  2017.388916015625,   2017.40283203125,
             2017.4166259765625, 2017.4305419921875, 2017.4444580078125,
             2017.4583740234375,   2017.47216796875,  2017.486083984375,
                         2017.5,  2017.513916015625,   2017.52783203125,
             2017.5416259765625, 2017.5555419921875, 2017.5694580078125,
             2017.5833740234375,   2017.59716796875,  2017.611083984375,
                       2017.625,  2017.638916015625,   2017.65283203125,
             2017.6666259765625, 2017.6805419921875, 2017.6944580078125,
             2017.7083740234375,   2017.72216796875,  2017.736083984375,
                        2017.75,  2017.763916015625,   2017.77783203125,
             2017.7916259765625, 2017.8055419921875, 2017.8194580078125,
             2017.8333740234375,   2017.84716796875,  2017.861083984375,
                       2017.875,  2017.888916015625,   2017.90283203125,
             2017.9166259765625, 2017.9305419921875, 2017.9444580078125,
             2017.9583740234375,   2017.97216796875,  2017.986083984375],
            dtype='float32', name='year'))
    • depth
      PandasIndex
      PandasIndex(Index([15.0], dtype='float32', name='depth'))
    • latitude
      PandasIndex
      PandasIndex(Index([              80.0,  79.66666666666667,  79.33333333333333,
                           79.0,  78.66666666666667,  78.33333333333333,
                           78.0,  77.66666666666667,  77.33333333333333,
                           77.0,
             ...
                          -77.0, -77.33333333333333, -77.66666666666667,
                          -78.0, -78.33333333333333, -78.66666666666667,
                          -79.0, -79.33333333333333, -79.66666666666667,
                          -80.0],
            dtype='float64', name='latitude', length=481))
    • longitude
      PandasIndex
      PandasIndex(Index([              20.0, 20.333333333333332, 20.666666666666664,
                           21.0, 21.333333333333332, 21.666666666666664,
                           22.0, 22.333333333333332, 22.666666666666664,
                           23.0,
             ...
                          417.0,  417.3333333333333,  417.6666666666667,
                          418.0,  418.3333333333333,  418.6666666666667,
                          419.0,  419.3333333333333,  419.6666666666667,
                          420.0],
            dtype='float64', name='longitude', length=1201))
  • VARIABLE :
    Ocean Surface Currents
    DATATYPE :
    1/72 YEAR Interval
    DATASUBTYPE :
    unfiltered
    GEORANGE :
    20 to 420 -80 to 80
    PERIOD :
    Jan.01,2017 to Dec.26,2017
    year :
    2017
    description :
    OSCAR Third Degree Sea Surface Velocity
    CREATION_DATE :
    03:24 05-Apr-2018
    version :
    2009.0
    source :
    Gary Lagerloef, ESR (lager@esr.org) and Kathleen Dohan, ESR (kdohan@esr.org)
    contact :
    Kathleen Dohan (kdohan@esr.org) or John T. Gunn (gunn@esr.org)
    company :
    Earth & Space Research, Seattle, WA
    reference :
    Bonjean F. and G.S.E. Lagerloef, 2002 ,Diagnostic model and analysis of the surface currents in the tropical Pacific ocean, J. Phys. Oceanogr., 32, 2,938-2,954
    note1 :
    Maximum Mask velocity is the geostrophic component at all points + any concurrent Ekman and buoyancy components
    note2 :
    Longitude extends from 20 E to 420 E to avoid a break in major ocean basins. Data repeats in overlap region.
In [18]:
# Select surface level (assuming depth=15.0 is the only level)
surface_currents2 = ds2.sel(depth=15.0)

# Calculate the annual mean for 2018
u_mean2 = surface_currents2['u'].mean(dim='time', skipna=True)
v_mean2 = surface_currents2['v'].mean(dim='time', skipna=True)

# Calculate magnitude of the currents
umag_g2 =(u_mean2**2 + v_mean2**2)**0.5

# Global plot
plt.figure(figsize=(14, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()  # Add coastlines
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')  # White background for land
ax.add_feature(cfeature.OCEAN, color='white')  # White background for oceans

quiver = plt.quiver(surface_currents2['longitude'], surface_currents2['latitude'], u_mean2, v_mean2, umag_g2, transform=ccrs.PlateCarree(), scale=25, pivot='middle', cmap='viridis')

ax.set_xticks(np.arange(-180, 181, 30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90, 91, 30), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x}°"))
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda y, _: f"{y}°"))
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")


plt.colorbar(quiver, label='Surface Current Magnitude (m/s)', orientation='vertical')
plt.title('Annual Mean Surface Currents (Global) - 2017')
plt.show()

# Regional subset for the northwestern Atlantic
region2 = surface_currents2.sel(latitude=slice(45, 25), longitude=slice(280, 300))
u_mean_region2 = region2['u'].mean(dim='time', skipna=True)
v_mean_region2 = region2['v'].mean(dim='time', skipna=True)
umag_r2 = (u_mean_region**2 + v_mean_region**2)**0.5
# Regional plot
plt.figure(figsize=(10, 6))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([280, 300, 25, 45], crs=ccrs.PlateCarree())
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='white')

quiver = plt.quiver(region2['longitude'], region2['latitude'], u_mean_region2, v_mean_region2, umag_r2, transform=ccrs.PlateCarree(), scale=10, pivot='middle', cmap='viridis')


cbar = plt.colorbar(quiver, label='Surface Current Magnitude (m/s)', orientation='vertical')


# Add gridlines with proper longitude and latitude labels
gl = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}

plt.title('Annual Mean Surface Currents (Northwestern Atlantic) - 2017')
plt.show()
No description has been provided for this image
No description has been provided for this image

Noticable change in the regional¶

#4. Plot the difference in annual means of the surface currents (both global and regional subset) between 2017 and 2018 as colored vectors. Describe the differences.¶

You can use plt.clim(vmin,vmax) to set your vmin and vmax outside of the quiver line, to make it easier to center your colorbar on 0. Be sure to pick a colormap that works well for this.¶
In [26]:
surface_currents_2017 = ds2.sel(depth=15.0)
surface_currents_2018 = ds.sel(depth=15.0)
u_mean_2017 = surface_currents_2017['u'].mean(dim='time', skipna=True)
v_mean_2017 = surface_currents_2017['v'].mean(dim='time', skipna=True)


u_diff = u_mean_2018 - u_mean_2017
v_diff = v_mean_2018 - v_mean_2017

# Compute the magnitude of the difference vectors
umag_diff = (u_diff**2 + v_diff**2)**0.5

max_diff = max(np.max(np.abs(u_diff)), np.max(np.abs(v_diff)))

# 1. Global plot
plt.figure(figsize=(14, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='white')

# Quiver plot for the global difference 
quiver = plt.quiver(surface_currents_2018['longitude'], surface_currents_2018['latitude'], u_diff, v_diff, umag_diff, transform=ccrs.PlateCarree(), scale=25, pivot='middle', cmap='coolwarm', clim=(-max_diff, max_diff))


ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")

plt.colorbar(quiver, label='Surface Current Difference Magnitude (m/s)', orientation='vertical')

# Gridlines with labels
gl = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}

# Title
plt.title('Difference in Annual Mean Surface Currents (Global) - 2018 vs 2017', fontsize=14)
plt.show()

# 2. Regional plot (Northwestern Atlantic)
plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([280, 300, 25, 45], crs=ccrs.PlateCarree())
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='white')

# Quiver plot for the regional difference 
quiver = plt.quiver(region_2018['longitude'], region_2018['latitude'], u_diff_region, v_diff_region, umag_diff_r, transform=ccrs.PlateCarree(), scale=10, pivot='middle', cmap='coolwarm', clim=(-max_diff, max_diff))

#Color bar
plt.colorbar(quiver, label='Surface Current Difference Magnitude (m/s)', orientation='vertical')

# Gridlines with labels
gl = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}

# Title
plt.title('Difference in Annual Mean Surface Currents (Northwestern Atlantic) - 2018 vs 2017', fontsize=14)
plt.show()
No description has been provided for this image
No description has been provided for this image

The above plots show the difference between 2018 currents and 2017 currents. The red is positive indicating stronger currents in 2018 comapered to 2017. In most cases curents stayed roughly the same. There is indication that the more prominent currents depicted in both years show a stronger average current in 2018. In the regional model there are red current vectors in in the top right of the plot showing a stronger gulf stream current during 2018.

#5. Plot the annual mean of the surface currents (both global and regional subset) for 2018 as colored streamlines.¶

In [35]:
subset_global = surface_currents_2018

# Meshgrid for longitude and latitude
lon_global, lat_global = np.meshgrid(subset_global.longitude - 360, np.flip(subset_global.latitude))

# Flips u, v for plotting
u_global = np.flipud(u_mean_2018.values)
v_global = np.flipud(v_mean_2018.values)
umag_global = np.sqrt(u_global**2 + v_global**2)

# Global streamline plot
plt.figure(figsize=(14, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()

# Coastlines and features
ax.coastlines(resolution='10m', linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='lightblue')

# Plot streamlines
strm = plt.streamplot(lon_global, lat_global, u_global, v_global, color=umag_global, linewidth=1, cmap='viridis', density=8, arrowsize=1, transform=ccrs.PlateCarree())

# Color bar
plt.colorbar(strm.lines, label='Surface Current Magnitude (m/s)')

ax.set_xticks(np.arange(-180, 181, 30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90, 91, 30), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x}°"))
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda y, _: f"{y}°"))
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")

# Title
plt.title('Global Annual Mean Surface Currents (2018)', fontsize=14)
plt.show()

# Regional streamline plot for Northwestern Atlantic
subset_regional = surface_currents_2018.sel(latitude=slice(45, 25), longitude=slice(280, 300))


lon_regional, lat_regional = np.meshgrid(subset_regional.longitude - 360, np.flip(subset_regional.latitude))

# Flip u, v for plotting
u_regional = np.flipud(subset_regional['u'].mean(dim='time', skipna=True).values)
v_regional = np.flipud(subset_regional['v'].mean(dim='time', skipna=True).values)
umag_regional = np.sqrt(u_regional**2 + v_regional**2)

# Regional streamline plot
plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([280 - 360, 300 - 360, 25, 45], crs=ccrs.PlateCarree())

# Coastlines and features
ax.coastlines(resolution='10m', linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, color='white')
ax.add_feature(cfeature.OCEAN, color='lightblue')

# Plot streamlines
strm = plt.streamplot(lon_regional, lat_regional, u_regional, v_regional, color=umag_regional, linewidth=1, cmap='viridis', density=8, arrowsize=1, transform=ccrs.PlateCarree())

# Color bar
plt.colorbar(strm.lines, label='Surface Current Magnitude (m/s)')

gl = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.5)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}

# Title
plt.title('Regional Annual Mean Surface Currents (Northwestern Atlantic, 2018)', fontsize=14)
plt.show()
No description has been provided for this image
No description has been provided for this image
In [ ]: